#ifndef _MATRIX_4X4_H_
#define _MATRIX_4X4_H_

/**
* @file matrix4x4.h
* July, 2010
* 
* 
* @author Antonio Miranda (amiranda@tecgraf.puc-rio.br )
* 
*/

#include <assert.h>

#include "point3d.h"

//
//////////////////////////////////////////////////////////////////////////
class matrix4x4
{
public:
  //matrix4x4(void);
  //~matrix4x4(void);

  inline matrix4x4() { setToIdentity(); }
  explicit matrix4x4(const double *values);
  inline matrix4x4(double m11, double m12, double m13, double m14,
                    double m21, double m22, double m23, double m24,
                    double m31, double m32, double m33, double m34,
                    double m41, double m42, double m43, double m44);

  matrix4x4(const double *values, int cols, int rows);

  inline void set (double m11, double m12, double m13, double m14,
                   double m21, double m22, double m23, double m24,
                   double m31, double m32, double m33, double m34,
                   double m41, double m42, double m43, double m44);

  inline const double& operator()(int row, int column) const;
  inline double& operator()(int row, int column);

  //inline QVector4D column(int index) const;
  //inline void setColumn(int index, const QVector4D& value);

  //inline QVector4D row(int index) const;
  //inline void setRow(int index, const QVector4D& value);

  inline bool isIdentity() const;
  inline void setToIdentity();

  inline void fill(double value);

  double    determinant() const;
  matrix4x4 inverted(bool *invertible = 0) const;
  matrix4x4 transposed() const;
  void      transposed (double m[4][4]) const;
  void      copy (double m[4][4]) const;
  //QMatrix3x3 normalMatrix() const;

  inline matrix4x4& operator+=(const matrix4x4& other);
  inline matrix4x4& operator-=(const matrix4x4& other);
  inline matrix4x4& operator*=(const matrix4x4& other);
  inline matrix4x4& operator*=(double factor);
  matrix4x4& operator/=(double divisor);
  inline bool operator==(const matrix4x4& other) const;
  inline bool operator!=(const matrix4x4& other) const;

  friend matrix4x4 operator+(const matrix4x4& m1, const matrix4x4& m2);
  friend matrix4x4 operator-(const matrix4x4& m1, const matrix4x4& m2);
  friend matrix4x4 operator*(const matrix4x4& m1, const matrix4x4& m2);

  friend Point3D operator*(const Point3D& point, const matrix4x4& matrix);
  friend matrix4x4 operator-(const matrix4x4& matrix);
  friend Point3D operator*(const matrix4x4& matrix, const Point3D& point);
  friend matrix4x4 operator*(double factor, const matrix4x4& matrix);
  friend matrix4x4 operator*(const matrix4x4& matrix, double factor);
  // friend matrix4x4 operator/(const matrix4x4 & matrix, double divisor);

  // friend inline bool qFuzzyCompare(const matrix4x4& m1, const matrix4x4& m2);

  void scale(double x, double y);
  void scale(double x, double y, double z);
  void scale(double factor);
  void translate(double x, double y);
  void translate(double x, double y, double z);
  void rotate(double angle, double x, double y, double z = 0.0f);

  //void ortho(const QRect& rect);
  //void ortho(const QRectF& rect);
  //void ortho(double left, double right, double bottom, double top, double nearPlane, double farPlane);
  //void frustum(double left, double right, double bottom, double top, double nearPlane, double farPlane);
  //void perspective(double angle, double aspect, double nearPlane, double farPlane);
  //void flipCoordinates();

  void copyDataTo(double *values) const;

  //QMatrix toAffine() const;
  //QTransform toTransform() const;
  //QTransform toTransform(double distanceToPlane) const;

  Point3D map(const Point3D& point) const;

  inline double *data();
  inline const double *data() const { return m[0]; }
  inline const double *constData() const { return m[0]; }

  void optimize();

  // create a base matrix where w is the vector from p1 to p2 
  static void baseMatrix(double p1[3], double p2[3], matrix4x4 &mat1, matrix4x4 &inv1); 


private:

  double m[4][4];          // Column-major order to match OpenGL.
  int flagBits;            // Flag bits from the enum below.

  enum 
  {
    Identity        = 0x0001,   // Identity matrix
    General         = 0x0002,   // General matrix, unknown contents
    Translation     = 0x0004,   // Contains a simple translation
    Scale           = 0x0008,   // Contains a simple scale
    Rotation        = 0x0010    // Contains a simple rotation
  };

  // Construct without initializing identity matrix.
  matrix4x4(int) { flagBits = General; }

  matrix4x4 orthonormalInverse() const;

  //void projectedRotate(double angle, double x, double y, double z);

};


inline matrix4x4::matrix4x4
(double m11, double m12, double m13, double m14,
 double m21, double m22, double m23, double m24,
 double m31, double m32, double m33, double m34,
 double m41, double m42, double m43, double m44)
{
  m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14;
  m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24;
  m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34;
  m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44;
  flagBits = General;
}

inline void matrix4x4::set
(double m11, double m12, double m13, double m14,
 double m21, double m22, double m23, double m24,
 double m31, double m32, double m33, double m34,
 double m41, double m42, double m43, double m44)
{
  m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14;
  m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24;
  m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34;
  m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44;
  flagBits = General;
}

inline const double& matrix4x4::operator()(int aRow, int aColumn) const
{
  assert(aRow >= 0 && aRow < 4 && aColumn >= 0 && aColumn < 4);
  return m[aColumn][aRow];
}

inline double& matrix4x4::operator()(int aRow, int aColumn)
{
  assert(aRow >= 0 && aRow < 4 && aColumn >= 0 && aColumn < 4);
  flagBits = General;
  return m[aColumn][aRow];
}

//inline QVector4D matrix4x4::column(int index) const
//{
//  assert(index >= 0 && index < 4);
//  return QVector4D(m[index][0], m[index][1], m[index][2], m[index][3]);
//}
//
//inline void matrix4x4::setColumn(int index, const QVector4D& value)
//{
//  assert(index >= 0 && index < 4);
//  m[index][0] = value.x();
//  m[index][1] = value.y();
//  m[index][2] = value.z();
//  m[index][3] = value.w();
//  flagBits = General;
//}
//
//inline QVector4D matrix4x4::row(int index) const
//{
//  assert(index >= 0 && index < 4);
//  return QVector4D(m[0][index], m[1][index], m[2][index], m[3][index]);
//}
//
//inline void matrix4x4::setRow(int index, const QVector4D& value)
//{
//  assert(index >= 0 && index < 4);
//  m[0][index] = value.x();
//  m[1][index] = value.y();
//  m[2][index] = value.z();
//  m[3][index] = value.w();
//  flagBits = General;
//}

// Q_GUI_EXPORT matrix4x4 operator/(const matrix4x4& matrix, double divisor);

inline bool matrix4x4::isIdentity() const
{
  if (flagBits == Identity)
    return true;
  if (m[0][0] != 1.0f || m[0][1] != 0.0f || m[0][2] != 0.0f)
    return false;
  if (m[0][3] != 0.0f || m[1][0] != 0.0f || m[1][1] != 1.0f)
    return false;
  if (m[1][2] != 0.0f || m[1][3] != 0.0f || m[2][0] != 0.0f)
    return false;
  if (m[2][1] != 0.0f || m[2][2] != 1.0f || m[2][3] != 0.0f)
    return false;
  if (m[3][0] != 0.0f || m[3][1] != 0.0f || m[3][2] != 0.0f)
    return false;
  return (m[3][3] == 1.0f);
}

inline void matrix4x4::setToIdentity()
{
  m[0][0] = 1.0f;
  m[0][1] = 0.0f;
  m[0][2] = 0.0f;
  m[0][3] = 0.0f;
  m[1][0] = 0.0f;
  m[1][1] = 1.0f;
  m[1][2] = 0.0f;
  m[1][3] = 0.0f;
  m[2][0] = 0.0f;
  m[2][1] = 0.0f;
  m[2][2] = 1.0f;
  m[2][3] = 0.0f;
  m[3][0] = 0.0f;
  m[3][1] = 0.0f;
  m[3][2] = 0.0f;
  m[3][3] = 1.0f;
  flagBits = Identity;
}

inline void matrix4x4::fill(double value)
{
  m[0][0] = value;
  m[0][1] = value;
  m[0][2] = value;
  m[0][3] = value;
  m[1][0] = value;
  m[1][1] = value;
  m[1][2] = value;
  m[1][3] = value;
  m[2][0] = value;
  m[2][1] = value;
  m[2][2] = value;
  m[2][3] = value;
  m[3][0] = value;
  m[3][1] = value;
  m[3][2] = value;
  m[3][3] = value;
  flagBits = General;
}

inline matrix4x4& matrix4x4::operator+=(const matrix4x4& other)
{
  m[0][0] += other.m[0][0];
  m[0][1] += other.m[0][1];
  m[0][2] += other.m[0][2];
  m[0][3] += other.m[0][3];
  m[1][0] += other.m[1][0];
  m[1][1] += other.m[1][1];
  m[1][2] += other.m[1][2];
  m[1][3] += other.m[1][3];
  m[2][0] += other.m[2][0];
  m[2][1] += other.m[2][1];
  m[2][2] += other.m[2][2];
  m[2][3] += other.m[2][3];
  m[3][0] += other.m[3][0];
  m[3][1] += other.m[3][1];
  m[3][2] += other.m[3][2];
  m[3][3] += other.m[3][3];
  flagBits = General;
  return *this;
}

inline matrix4x4& matrix4x4::operator-=(const matrix4x4& other)
{
  m[0][0] -= other.m[0][0];
  m[0][1] -= other.m[0][1];
  m[0][2] -= other.m[0][2];
  m[0][3] -= other.m[0][3];
  m[1][0] -= other.m[1][0];
  m[1][1] -= other.m[1][1];
  m[1][2] -= other.m[1][2];
  m[1][3] -= other.m[1][3];
  m[2][0] -= other.m[2][0];
  m[2][1] -= other.m[2][1];
  m[2][2] -= other.m[2][2];
  m[2][3] -= other.m[2][3];
  m[3][0] -= other.m[3][0];
  m[3][1] -= other.m[3][1];
  m[3][2] -= other.m[3][2];
  m[3][3] -= other.m[3][3];
  flagBits = General;
  return *this;
}

inline matrix4x4& matrix4x4::operator*=(const matrix4x4& other)
{
  if (flagBits == Identity) {
    *this = other;
    return *this;
  } else if (other.flagBits == Identity) {
    return *this;
  } else {
    *this = *this * other;
    return *this;
  }
}

inline matrix4x4& matrix4x4::operator*=(double factor)
{
  m[0][0] *= factor;
  m[0][1] *= factor;
  m[0][2] *= factor;
  m[0][3] *= factor;
  m[1][0] *= factor;
  m[1][1] *= factor;
  m[1][2] *= factor;
  m[1][3] *= factor;
  m[2][0] *= factor;
  m[2][1] *= factor;
  m[2][2] *= factor;
  m[2][3] *= factor;
  m[3][0] *= factor;
  m[3][1] *= factor;
  m[3][2] *= factor;
  m[3][3] *= factor;
  flagBits = General;
  return *this;
}

inline bool matrix4x4::operator==(const matrix4x4& other) const
{
  return m[0][0] == other.m[0][0] &&
    m[0][1] == other.m[0][1] &&
    m[0][2] == other.m[0][2] &&
    m[0][3] == other.m[0][3] &&
    m[1][0] == other.m[1][0] &&
    m[1][1] == other.m[1][1] &&
    m[1][2] == other.m[1][2] &&
    m[1][3] == other.m[1][3] &&
    m[2][0] == other.m[2][0] &&
    m[2][1] == other.m[2][1] &&
    m[2][2] == other.m[2][2] &&
    m[2][3] == other.m[2][3] &&
    m[3][0] == other.m[3][0] &&
    m[3][1] == other.m[3][1] &&
    m[3][2] == other.m[3][2] &&
    m[3][3] == other.m[3][3];
}

inline bool matrix4x4::operator!=(const matrix4x4& other) const
{
  return m[0][0] != other.m[0][0] ||
    m[0][1] != other.m[0][1] ||
    m[0][2] != other.m[0][2] ||
    m[0][3] != other.m[0][3] ||
    m[1][0] != other.m[1][0] ||
    m[1][1] != other.m[1][1] ||
    m[1][2] != other.m[1][2] ||
    m[1][3] != other.m[1][3] ||
    m[2][0] != other.m[2][0] ||
    m[2][1] != other.m[2][1] ||
    m[2][2] != other.m[2][2] ||
    m[2][3] != other.m[2][3] ||
    m[3][0] != other.m[3][0] ||
    m[3][1] != other.m[3][1] ||
    m[3][2] != other.m[3][2] ||
    m[3][3] != other.m[3][3];
}

inline matrix4x4 operator+(const matrix4x4& m1, const matrix4x4& m2)
{
  matrix4x4 m(1);
  m.m[0][0] = m1.m[0][0] + m2.m[0][0];
  m.m[0][1] = m1.m[0][1] + m2.m[0][1];
  m.m[0][2] = m1.m[0][2] + m2.m[0][2];
  m.m[0][3] = m1.m[0][3] + m2.m[0][3];
  m.m[1][0] = m1.m[1][0] + m2.m[1][0];
  m.m[1][1] = m1.m[1][1] + m2.m[1][1];
  m.m[1][2] = m1.m[1][2] + m2.m[1][2];
  m.m[1][3] = m1.m[1][3] + m2.m[1][3];
  m.m[2][0] = m1.m[2][0] + m2.m[2][0];
  m.m[2][1] = m1.m[2][1] + m2.m[2][1];
  m.m[2][2] = m1.m[2][2] + m2.m[2][2];
  m.m[2][3] = m1.m[2][3] + m2.m[2][3];
  m.m[3][0] = m1.m[3][0] + m2.m[3][0];
  m.m[3][1] = m1.m[3][1] + m2.m[3][1];
  m.m[3][2] = m1.m[3][2] + m2.m[3][2];
  m.m[3][3] = m1.m[3][3] + m2.m[3][3];
  return m;
}

inline matrix4x4 operator-(const matrix4x4& m1, const matrix4x4& m2)
{
  matrix4x4 m(1);
  m.m[0][0] = m1.m[0][0] - m2.m[0][0];
  m.m[0][1] = m1.m[0][1] - m2.m[0][1];
  m.m[0][2] = m1.m[0][2] - m2.m[0][2];
  m.m[0][3] = m1.m[0][3] - m2.m[0][3];
  m.m[1][0] = m1.m[1][0] - m2.m[1][0];
  m.m[1][1] = m1.m[1][1] - m2.m[1][1];
  m.m[1][2] = m1.m[1][2] - m2.m[1][2];
  m.m[1][3] = m1.m[1][3] - m2.m[1][3];
  m.m[2][0] = m1.m[2][0] - m2.m[2][0];
  m.m[2][1] = m1.m[2][1] - m2.m[2][1];
  m.m[2][2] = m1.m[2][2] - m2.m[2][2];
  m.m[2][3] = m1.m[2][3] - m2.m[2][3];
  m.m[3][0] = m1.m[3][0] - m2.m[3][0];
  m.m[3][1] = m1.m[3][1] - m2.m[3][1];
  m.m[3][2] = m1.m[3][2] - m2.m[3][2];
  m.m[3][3] = m1.m[3][3] - m2.m[3][3];
  return m;
}

inline matrix4x4 operator*(const matrix4x4& m1, const matrix4x4& m2)
{
  if (m1.flagBits == matrix4x4::Identity)
    return m2;
  else if (m2.flagBits == matrix4x4::Identity)
    return m1;

  matrix4x4 m(1);
  double sum = 0.0;
  for ( int c = 0 ; c < 4 ; c++ )
  {
    for ( int d = 0 ; d < 4 ; d++ )
    {
      for ( int k = 0 ; k < 4 ; k++ )
      {
        sum = sum + m1.m[c][k]*m2.m[k][d];
      }

      m.m[c][d] = sum;
      sum = 0.0;
    }
  }

  return m;
}

inline Point3D operator*(const Point3D& vector, const matrix4x4& matrix)
{
  double x, y, z, w;
  x = vector.getX() * matrix.m[0][0] +
    vector.getY() * matrix.m[0][1] +
    vector.getZ() * matrix.m[0][2] +
    matrix.m[0][3];
  y = vector.getX() * matrix.m[1][0] +
    vector.getY() * matrix.m[1][1] +
    vector.getZ() * matrix.m[1][2] +
    matrix.m[1][3];
  z = vector.getX() * matrix.m[2][0] +
    vector.getY() * matrix.m[2][1] +
    vector.getZ() * matrix.m[2][2] +
    matrix.m[2][3];
  w = vector.getX() * matrix.m[3][0] +
    vector.getY() * matrix.m[3][1] +
    vector.getZ() * matrix.m[3][2] +
    matrix.m[3][3];
  //if (w == 1.0f)
  //  return Point3D(x, y, z);
  //else
  //  return Point3D(x / w, y / w, z / w);
}

inline Point3D operator*(const matrix4x4& matrix, const Point3D& vector)
{
  double x, y, z;

  if (matrix.flagBits == matrix4x4::Identity) 
  {
    return vector;
  } 
  else if (matrix.flagBits == matrix4x4::Translation) 
  {
    return Point3D(vector.getX() + matrix.m[3][0],
      vector.getY() + matrix.m[3][1],
      vector.getZ() + matrix.m[3][2]);
  } 
  else if (matrix.flagBits == (matrix4x4::Translation | matrix4x4::Scale)) 
  {
    return Point3D(vector.getX() * matrix.m[0][0] + matrix.m[3][0],
      vector.getY() * matrix.m[1][1] + matrix.m[3][1],
      vector.getZ() * matrix.m[2][2] + matrix.m[3][2]);
  } 
  else if (matrix.flagBits == matrix4x4::Scale) 
  {
    return Point3D(vector.getX() * matrix.m[0][0],
      vector.getY() * matrix.m[1][1],
      vector.getZ() * matrix.m[2][2]);
  } 
  else 
  {
    x = vector.getX() * matrix.m[0][0] +
      vector.getY() * matrix.m[0][1] +
      vector.getZ() * matrix.m[0][2] +
      matrix.m[0][3];
    y = vector.getX() * matrix.m[1][0] +
      vector.getY() * matrix.m[1][1] +
      vector.getZ() * matrix.m[1][2] +
      matrix.m[1][3];
    z = vector.getX() * matrix.m[2][0] +
      vector.getY() * matrix.m[2][1] +
      vector.getZ() * matrix.m[2][2] +
      matrix.m[2][3];
    //w = vector.getX() * matrix.m[0][3] +
    //  vector.getY() * matrix.m[1][3] +
    //  vector.getZ() * matrix.m[2][3] +
    //  matrix.m[3][3];
    return Point3D(x, y, z);
    //if (w == 1.0f)
    //  return Point3D(x, y, z);
    //else
    //  return Point3D(x / w, y / w, z / w);
  }
}





inline matrix4x4 operator-(const matrix4x4& matrix)
{
  matrix4x4 m(1);
  m.m[0][0] = -matrix.m[0][0];
  m.m[0][1] = -matrix.m[0][1];
  m.m[0][2] = -matrix.m[0][2];
  m.m[0][3] = -matrix.m[0][3];
  m.m[1][0] = -matrix.m[1][0];
  m.m[1][1] = -matrix.m[1][1];
  m.m[1][2] = -matrix.m[1][2];
  m.m[1][3] = -matrix.m[1][3];
  m.m[2][0] = -matrix.m[2][0];
  m.m[2][1] = -matrix.m[2][1];
  m.m[2][2] = -matrix.m[2][2];
  m.m[2][3] = -matrix.m[2][3];
  m.m[3][0] = -matrix.m[3][0];
  m.m[3][1] = -matrix.m[3][1];
  m.m[3][2] = -matrix.m[3][2];
  m.m[3][3] = -matrix.m[3][3];
  return m;
}

inline matrix4x4 operator*(double factor, const matrix4x4& matrix)
{
  matrix4x4 m(1);
  m.m[0][0] = matrix.m[0][0] * factor;
  m.m[0][1] = matrix.m[0][1] * factor;
  m.m[0][2] = matrix.m[0][2] * factor;
  m.m[0][3] = matrix.m[0][3] * factor;
  m.m[1][0] = matrix.m[1][0] * factor;
  m.m[1][1] = matrix.m[1][1] * factor;
  m.m[1][2] = matrix.m[1][2] * factor;
  m.m[1][3] = matrix.m[1][3] * factor;
  m.m[2][0] = matrix.m[2][0] * factor;
  m.m[2][1] = matrix.m[2][1] * factor;
  m.m[2][2] = matrix.m[2][2] * factor;
  m.m[2][3] = matrix.m[2][3] * factor;
  m.m[3][0] = matrix.m[3][0] * factor;
  m.m[3][1] = matrix.m[3][1] * factor;
  m.m[3][2] = matrix.m[3][2] * factor;
  m.m[3][3] = matrix.m[3][3] * factor;
  return m;
}

inline matrix4x4 operator*(const matrix4x4& matrix, double factor)
{
  matrix4x4 m(1);
  m.m[0][0] = matrix.m[0][0] * factor;
  m.m[0][1] = matrix.m[0][1] * factor;
  m.m[0][2] = matrix.m[0][2] * factor;
  m.m[0][3] = matrix.m[0][3] * factor;
  m.m[1][0] = matrix.m[1][0] * factor;
  m.m[1][1] = matrix.m[1][1] * factor;
  m.m[1][2] = matrix.m[1][2] * factor;
  m.m[1][3] = matrix.m[1][3] * factor;
  m.m[2][0] = matrix.m[2][0] * factor;
  m.m[2][1] = matrix.m[2][1] * factor;
  m.m[2][2] = matrix.m[2][2] * factor;
  m.m[2][3] = matrix.m[2][3] * factor;
  m.m[3][0] = matrix.m[3][0] * factor;
  m.m[3][1] = matrix.m[3][1] * factor;
  m.m[3][2] = matrix.m[3][2] * factor;
  m.m[3][3] = matrix.m[3][3] * factor;
  return m;
}

//inline bool qFuzzyCompare(const matrix4x4& m1, const matrix4x4& m2)
//{
//  return qFuzzyCompare(m1.m[0][0], m2.m[0][0]) &&
//    qFuzzyCompare(m1.m[0][1], m2.m[0][1]) &&
//    qFuzzyCompare(m1.m[0][2], m2.m[0][2]) &&
//    qFuzzyCompare(m1.m[0][3], m2.m[0][3]) &&
//    qFuzzyCompare(m1.m[1][0], m2.m[1][0]) &&
//    qFuzzyCompare(m1.m[1][1], m2.m[1][1]) &&
//    qFuzzyCompare(m1.m[1][2], m2.m[1][2]) &&
//    qFuzzyCompare(m1.m[1][3], m2.m[1][3]) &&
//    qFuzzyCompare(m1.m[2][0], m2.m[2][0]) &&
//    qFuzzyCompare(m1.m[2][1], m2.m[2][1]) &&
//    qFuzzyCompare(m1.m[2][2], m2.m[2][2]) &&
//    qFuzzyCompare(m1.m[2][3], m2.m[2][3]) &&
//    qFuzzyCompare(m1.m[3][0], m2.m[3][0]) &&
//    qFuzzyCompare(m1.m[3][1], m2.m[3][1]) &&
//    qFuzzyCompare(m1.m[3][2], m2.m[3][2]) &&
//    qFuzzyCompare(m1.m[3][3], m2.m[3][3]);
//}


inline Point3D matrix4x4::map(const Point3D& point) const
{
  return *this * point;
}

//inline Point3D matrix4x4::map(const Point3D& vector) const
//{
//  if (flagBits == Identity || flagBits == Translation) {
//    return vector;
//  } else if (flagBits == Scale || flagBits == (Translation | Scale)) {
//    return Point3D(vector.getX() * m[0][0],
//      vector.getY() * m[1][1],
//      vector.getZ() * m[2][2]);
//  } else {
//    return Point3D(vector.getX() * m[0][0] +
//      vector.getY() * m[1][0] +
//      vector.getZ() * m[2][0],
//      vector.getX() * m[0][1] +
//      vector.getY() * m[1][1] +
//      vector.getZ() * m[2][1],
//      vector.getX() * m[0][2] +
//      vector.getY() * m[1][2] +
//      vector.getZ() * m[2][2]);
//  }
//}


inline double *matrix4x4::data()
{
  // We have to assume that the caller will modify the matrix elements,
  // so we flip it over to "General" mode.
  flagBits = General;
  return m[0];
}

#endif